
% HHVR_EM1a_2_types.m is linked to HHVR_EM1b_2_types. 
 
function  [parameters_EM   z Z  LOGLIK loglik   PI EM_increasing] = HHVR_EM1a_2_typesBOOTSTRAP(initial_guess,who,DATA_ALL )
iterations=333;
convergence = 2 ;   

global DATA choices y y_short X M zz individuals valid_rounds
 
LogLik = [2222 1 ] ;
LOGLIK=[]; 
bbeta = initial_guess  ;
 
pi_u = 1/2; pi_e = 1/2; 
s=size(DATA,1)/choices;     % individuals x choices
nume1=zeros(s,1);nume2=zeros(s,1);nume3=zeros(s,1) ;nume4=zeros(s,1) ;nume5=zeros(s,1) ; 
numu1=zeros(s,1);numu2=zeros(s,1);numu3=zeros(s,1) ;numu4=zeros(s,1) ;numu5=zeros(s,1) ; 
prob_e1=zeros(s,1);prob_e2=zeros(s,1);prob_e3=zeros(s,1); prob_e4=zeros(s,1);prob_e5=zeros(s,1); 
prob_u1=zeros(s,1);prob_u2=zeros(s,1);prob_u3=zeros(s,1); prob_u4=zeros(s,1);prob_u5=zeros(s,1);  

denu= zeros(s,1); dene= zeros(s,1);    denu_x1= zeros(s,1); dene_x1= zeros(s,1);  

iter=0;
s=size(DATA,1)/choices;     % individuals x choices
numc1=zeros(s,1);numc2=zeros(s,1);numc3=zeros(s,1) ;numc4=zeros(s,1) ;numc5=zeros(s,1) ;numc6=zeros(s,1) ;
nump1=zeros(s,1);nump2=zeros(s,1);nump3=zeros(s,1); nump4=zeros(s,1);nump5=zeros(s,1);nump6=zeros(s,1); 

prob_c1=zeros(s,1);prob_c2=zeros(s,1);prob_c3=zeros(s,1); prob_c4=zeros(s,1);prob_c5=zeros(s,1);prob_c6=zeros(s,1); 
prob_p1=zeros(s,1);prob_p2=zeros(s,1);prob_p3=zeros(s,1); prob_p4=zeros(s,1);prob_p5=zeros(s,1); prob_p6=zeros(s,1);  
denc= zeros(s,1); denp= zeros(s,1); den_cx1= zeros(s,1);denp_x1= zeros(s,1);

nump7=zeros(s,1); 
prob_c7=zeros(s,1); numc7=zeros(s,1) ;  prob_p7=zeros(s,1);prob_p8=zeros(s,1);
Z=[];
for t=1:iterations
            while abs(LogLik(end-1)-LogLik(end)) > convergence
                   DIFF_LOG_LIK= (LogLik(end-1)-LogLik(end) ) 
                   DENE=[]; NUME=[];
                    
    if t <=2
        EM_increasing=0;
    end
    while t > 3
    if (LogLik(end-2)-LogLik(end)>0 ) 
    EM_increasing = who;
    elseif (LogLik(end-2)-LogLik(end)<=0 ) 
        EM_increasing=EM_increasing;
    end
    end
    
for i=1:(size(DATA,1))/choices % ie from 1 to number of observations per individual (48) 
 
    for h = 1:5 
         if  (X(i*choices-h+1,1)*bbeta(5)) > 707
         X(i*choices-h+1,1)  = 707/bbeta(5);
         %disp 'oops' 
         %pause
         end
    end
 
     nume1(i) =  exp(X(i*choices-4,1)*bbeta(5)    ) ;
     nume2(i) =  exp(X(i*choices-3,1)*bbeta(5)  + 1*bbeta(1)  ) ;
     nume3(i) =  exp(X(i*choices-2,1)*bbeta(5)  + 1*bbeta(2)  ) ; 
 
     nume4(i) =  exp(X(i*choices-1,1)*bbeta(5)  + 1*bbeta(3)  ) ;
     nume5(i) =  exp(X(i*choices  ,1)*bbeta(5)  + 1*bbeta(4)  ) ;
    
     dene(i)   =  nume1 (i)  +  nume2 (i) + nume3 (i)+ nume4 (i)+ nume5 (i);
     dene_x1(i)   =  nume1(i)*X(i*choices-4,1)  +  nume2(i)*X(i*choices-3,1) + nume3(i)*X(i*choices-2,1)+ ... 
                     nume4(i)*X(i*choices-1,1) +   nume5(i)*X(i*choices  ,1 ) ;
     prob_e1(i)     =  nume1 (i) /  dene(i)  ;     prob_e2(i)     =  nume2 (i) /  dene(i)  ;
     prob_e3(i)     =  nume3 (i) /  dene(i)  ;     prob_e4(i)     =  nume4 (i) /  dene(i)  ;
     prob_e5(i)     =  nume5 (i) /  dene(i)  ;  
    
     numu1 (i) =  exp(X(i*choices-4,2)*bbeta(6)    ) ;
     numu2 (i) =  exp(X(i*choices-3,2)*bbeta(6)  + 1*bbeta(1)  ) ;
     numu3 (i) =  exp(X(i*choices-2,2)*bbeta(6)  + 1*bbeta(2)  ) ;
     numu4 (i) =  exp(X(i*choices-1,2)*bbeta(6)  + 1*bbeta(3)  ) ;
     numu5 (i) =  exp(X(i*choices  ,2)*bbeta(6)  + 1*bbeta(4)  ) ;
     denu(i)   =  numu1 (i)  +  numu2 (i) + numu3 (i)+ numu4 (i)+ numu5 (i) ;
     denu_x1(i)   =  numu1(i)*X(i*choices-4,2)  +  numu2(i)*X(i*choices-3,2) + numu3(i)*X(i*choices-2,2)+ ... 
                     numu4(i)*X(i*choices-1,2) +   numu5(i)*X(i*choices  ,2 ) ;
     prob_u1(i)     =  numu1 (i) /  denu(i)  ;     prob_u2(i)     =  numu2 (i) /  denu(i)  ;
     prob_u3(i)     =  numu3 (i) /  denu(i)  ;     prob_u4(i)     =  numu4 (i) /  denu(i)  ;
     prob_u5(i)     =  numu5 (i) /  denu(i)  ;  
     DENE = [DENE ; dene(i)];
     NUME = [NUME ; nume1(i) nume2(i) nume3(i) nume4(i) nume5(i)];
end
 
densitye=zeros(s,1)'; densityu=zeros(s,1)';  
        
 
for i=1:(size(DATA,1))/choices
    if y_short(i)==1
        densitye(i)= prob_e1(i)  ;       densityu(i)= prob_u1(i) ;
    elseif y_short(i)==2
        densitye(i)= prob_e2(i) ;        densityu(i)= prob_u2(i) ; 
    elseif y_short(i)==3
        densitye(i)= prob_e3(i) ;        densityu(i)= prob_u3(i) ; 
    elseif y_short(i)==4
        densitye(i)= prob_e4(i) ;        densityu(i)= prob_u4(i) ; 
    elseif y_short(i)==5
        densitye(i)= prob_e5(i) ;        densityu(i)= prob_u5(i) ; 
     end
end

density=[densitye ; densityu    ]   ;        % 3  x individuals and they are all different
ze=zeros((size(DATA,1))/choices,1);                % individuals  x  1
zu=zeros((size(DATA,1))/choices,1);                % individuals  x  1 

for i=1:(size(DATA,1))/choices
    ze(i) = ( pi_e*densitye(i)) / (pi_e*densitye(i) + pi_u*densityu(i)    );
    zu(i) = ( pi_u*densityu(i))    / (pi_e*densitye(i) + pi_u*densityu(i)   ); 
 end
Z=[Z ze zu  ];
   
            loglik=[];Loglikvuong = [];
            for i=1:s
                loglik(i) =  ze(i)*log(densitye(i)) + zu(i)*log(densityu(i))    ;
                Loglikvuong = [Loglikvuong ; loglik(i)];
            end
            Loglik = sum(loglik);
            LogLik = [LogLik Loglik]  ;
            DIFF_LOG_LIK= (LogLik(end-1)-LogLik(end) );
            pi_e = mean(ze); pi_u = mean(zu);  
            LOGLIK = [LOGLIK ; LogLik(end)  DIFF_LOG_LIK];
            LogLik(end);
            PI = [pi_e pi_u  ] 
            zz = [ze zu ] ; 
            inputs = [bbeta(1:6) ]  ;
            if abs(inputs (5)) > 700
                inputs(5)=700
            end
 
[bbeta ] = fsolve(@(bbeta)    HHVR_EM1b_2_types(bbeta), inputs,optimset('maxiter',20000       , 'MaxFunEvals',40000     ) ) ; % not [bbeta,fval,a,b]


bbeta(1:4);
CONSTANT =  [bbeta(1) bbeta(2)  bbeta(3) bbeta(4)    ];
if  bbeta(5)>700
    bbeta(5)=700;
end
BBETA = [bbeta(5) bbeta(6)   ]  ;
 
iter=iter+1;

            end %while
end

LOGLIK=LOGLIK' 
LogLik =  LogLik';
z=zz;
 
if size(DATA,2)>12
id_agents = reshape(DATA(:,1)',5,numel(DATA(:,1))/5) ;
id_agents = id_agents(1,:);
disjoint = reshape(DATA(:,13)',5,numel(DATA(:,13))/5) ;
disjoint = disjoint(1,:);
size(disjoint)
LogLik =  LogLik';
z=[ zz id_agents' disjoint'];
else
    z = zz;
end
 
total_vote_type = [0 0 0 ];
for i = 1: size(DATA,1)
    if DATA(i,2)== 1 && DATA(i,3) ==1
        total_vote_type(1) = total_vote_type(1) + 1;
    elseif DATA(i,2)== 1 && DATA(i,4) ==1
        total_vote_type(2) = total_vote_type(2) + 1;
    elseif DATA(i,2)== 1 && DATA(i,5) ==1
        total_vote_type(3) = total_vote_type(3) + 1;
    end
end
parameters_EM = [bbeta ];



 



